{
 "metadata": {
  "name": "",
  "signature": "sha256:782a2ac76756b16b8c29e6adba08c9019b65be21fff8cec4257b5e0204119fb5"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Blind Source Separation on Images with Shogun"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "by Kevin Hughes"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook illustrates <a href=\"http://en.wikipedia.org/wiki/Blind_signal_separation\">Blind Source Seperation</a>(BSS) on images using <a href=\"http://en.wikipedia.org/wiki/Independent_component_analysis\">Independent Component Analysis</a> (ICA) in Shogun. This is very similar to the <a href=\"http://www.shogun-toolbox.org/static/notebook/current/bss_audio.html\">BSS audio notebook</a> except that here we have used images instead of audio signals."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The first step is to load 2 images from the Shogun data repository:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# change to the shogun-data directory\n",
      "import os\n",
      "import os\nSHOGUN_DATA_DIR=os.getenv('SHOGUN_DATA_DIR', '../../../data')\n",
      "os.chdir(os.path.join(SHOGUN_DATA_DIR, 'ica'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import Image\n",
      "import numpy as np\n",
      "\n",
      "# Load Images as grayscale images and convert to numpy arrays\n",
      "s1 = np.asarray(Image.open(\"lena.jpg\").convert('L'))\n",
      "s2 = np.asarray(Image.open(\"monalisa.jpg\").convert('L'))\n",
      "\n",
      "# Save Image Dimensions\n",
      "# we'll need these later for reshaping the images\n",
      "rows = s1.shape[0]\n",
      "cols = s1.shape[1]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Displaying the images using pylab:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import pylab as pl\n",
      "\n",
      "# Show Images\n",
      "f,(ax1,ax2) = pl.subplots(1,2)\n",
      "ax1.imshow(s1, cmap=pl.gray()) # set the color map to gray, only needs to be done once!\n",
      "ax2.imshow(s2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In our previous ICA examples the input data or source signals were already 1D but these images are obviously 2D. One common way to handle this case is to simply \"flatten\" the 2D image matrix into a 1D row vector. The same idea can also be applied to 3D data, for example a 3 channel RGB image can be converted a row vector by reshaping each 2D channel into a row vector and then placing them after each other length wise.\n",
      "\n",
      "Lets prep the data:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Convert Images to row vectors\n",
      "# and stack into a Data Matrix\n",
      "S = np.c_[s1.flatten(), s2.flatten()].T"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It is pretty easy using a nice library like numpy.\n",
      "\n",
      "Next we need to mix our source signals together. We do this exactly the same way we handled the audio data - take a look!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Mixing Matrix\n",
      "A = np.array([[1, 0.5], [0.5, 1]])\n",
      "\n",
      "# Mix Signals\n",
      "X = np.dot(A,S)\n",
      "\n",
      "# Show Images\n",
      "f,(ax1,ax2) = pl.subplots(1,2)\n",
      "ax1.imshow(X[0,:].reshape(rows,cols))\n",
      "ax2.imshow(X[1,:].reshape(rows,cols))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Notice how we had to reshape from a 1D row vector back into a 2D matrix of the correct shape. There is also another nuance that I would like to mention here: pylab is actually doing quite a lot for us here that you might not be aware of. It does a pretty good job determining the value range of the image to be shown and then it applies the color map. Many other libraries (for example OpenCV's highgui) won't be this helpful and you'll need to remember to scale the image appropriately on your own before trying to display it. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now onto the exciting step, unmixing the images using ICA! Again this step is the same as when using Audio data. Again we need to reshape the images before viewing them and an additional nuance was to add the *-1 to the first separated signal. I did this after viewing the result the first time as the image was clearly inversed, this can happen because ICA can't necessarily capture the correct phase."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from shogun  import RealFeatures\n",
      "from shogun import Jade\n",
      "\n",
      "mixed_signals = RealFeatures(X)\n",
      "\n",
      "# Separating\n",
      "jade = Jade()\n",
      "signals = jade.apply(mixed_signals)\n",
      "S_ = signals.get_feature_matrix()\n",
      "\n",
      "# Show Images\n",
      "f,(ax1,ax2) = pl.subplots(1,2)\n",
      "ax1.imshow(S_[0,:].reshape(rows,cols) *-1)\n",
      "ax2.imshow(S_[1,:].reshape(rows,cols))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And that's all there is to it!"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}
